#include "core/full-dyn.hpp"
#include "core/cell-dyn-mt.hpp"
#include "core/verlet.hpp"
#include "core/timings.hpp"
#include "core/setup.hpp"
#include "core/ljc.hpp"
#include "core/fcc.hpp"
#include "core/constants.hpp"

#include "yocto/auto-ptr.hpp"
#include "yocto/lua/lua-state.hpp"
#include "yocto/lua/lua-config.hpp"
#include "yocto/exceptions.hpp"
#include "yocto/string/vfs-utils.hpp"
#include "yocto/eta.hpp"
#include "yocto/duration.hpp"

#include <cerrno>
#include <iomanip>

using namespace FastMD;

//==============================================================================
// format time
//==============================================================================
static inline const string time_fmt( const double ns )
{
	const  duration d(ns);
	return vformat("%02ud:%02u:%02u:%04.1fs", d.d, d.h, d.m, d.s );
}

//==============================================================================
// progress callback
//==============================================================================
static inline void progress( const eta &prg, const double nops )
{
	const string s_done = time_fmt( prg.time_done );
	const string s_left = time_fmt( prg.time_left );
	fprintf( stderr, "%4.2f%% in %s | ETA= %s | @%5.1f MOPS   \r", 
			prg.ratio_done*100.0, s_done.c_str(), s_left.c_str(), 1e-6*nops );
	fflush( stderr );
}

//==============================================================================
// display information
//==============================================================================
#define __SHOW_WIDTH 20
static inline size_t __width( const char *txt ) throw()
{
	assert(txt);
	const size_t n = strlen( txt );
	if( n >= __SHOW_WIDTH ) return 0;
	return __SHOW_WIDTH-n;
}
#define __SHOW(PARAM) std::cerr << "-- " <<  #PARAM << std::setw( __width(#PARAM) ) << " = " << PARAM << std::endl



int main(int argc, char *argv[])
{
	const char *program = _vfs::get_base_name( argv[0] );
	try
	{
		if( argc < 2 )
			throw exception("Usage: %s config.lua [lua code...]", program);
		const string config_name = argv[1];
		const string config_path = _vfs::get_file_dir( config_name );
		__SHOW(config_name);
		__SHOW(config_path);
		
		////////////////////////////////////////////////////////////////////////
		//
		// Resources to be Allocated
		//
		////////////////////////////////////////////////////////////////////////
		auto_ptr<Step>      step;
		auto_ptr<Frame>     frame;
		auto_ptr<Potential> potential; 
		auto_ptr<Dynamics>  dynamics;
		auto_ptr<Timings>   timings;
		auto_ptr<string>    pOutput;
		Real                t_ini = 0;
		Real                t_sim = 0;
		bool                isothermal = false;
		SRand();
		
		////////////////////////////////////////////////////////////////////////
		//
		// Read parameters
		//
		////////////////////////////////////////////////////////////////////////
		std::cerr << "-------- Parameters --------" << std::endl;
		{
			Lua::State VM;
			lua_State *L = VM();
			Lua::Config::DoFile( L, config_name );
			for( int i=2; i < argc; ++i )
			{
				const string cmd = argv[i];
				Lua::Config::DoString( L, cmd );
			}
			//==================================================================
			//
			// Read units
			//
			//==================================================================
			const string units = Setup::GetString(L,"units");
			while( true )
			{
				if( units == "reduced" )
				{
					Constants::Reduced();
					break;
				}
				
				if( units == "physical" )
				{
					Constants::Physical();
					break;
				}
				
				throw exception("Unknown units %s", units.c_str());
			}
			
			//==================================================================
			//
			// Read timings
			//
			//==================================================================
			const Real dt_run  = Setup::GetReal(L,"dt_run");  if( dt_run  <= 0 ) throw libc::exception( EINVAL, "dt_run<=0" );  
			const Real dt_save = Setup::GetReal(L,"dt_save"); if( dt_save <= 0 ) throw libc::exception( EINVAL, "dt_save<=0" );
			timings.reset( new Timings(dt_run,dt_save) );
			__SHOW(timings->dt_run);
			__SHOW(timings->dt_save);
			__SHOW(timings->every);
			
			//==================================================================
			//
			// Prepare Step
			//
			//==================================================================
			const string step_name = Setup::GetString(L,"step");
			while( true )
			{
				if( step_name == "verlet" )
				{
					step.reset( new Verlet( timings->dt_run ) );
					break;
				}
				
				throw libc::exception( EINVAL, "unknowm step '%s'", step_name.c_str() );
			}
			assert(step.__get() != NULL);
			
			//==================================================================
			//
			// Prepare Initial Frame
			//
			//==================================================================
			const string frame_type = Setup::GetString(L,"frame_type");
			__SHOW(frame_type);
			while(true)
			{
				
				if( frame_type == "generate" )
				{
					
					//----------------------------------------------------------
					//-- get the output name
					//----------------------------------------------------------
					const string output = Setup::GetString( L, "output" );
					if( length_of( _vfs::get_base_name(output) ) <= 0 )
						throw libc::exception( EINVAL,  "invalide output '%s'", output.c_str() );
					const string output_name = config_path + output;
					std::cerr << "saving in '" << output_name << "'" << std::endl;
					pOutput.reset( new string( output_name ) );
					
					
					//----------------------------------------------------------
					//-- get the box
					//----------------------------------------------------------
					const V3D gen_box = Setup::GetV3D( L, "gen_box" ); 
					if( gen_box.x <= 0 || gen_box.y <=0 || gen_box.z <= 0 ) throw libc::exception( EINVAL, "invalid gen_box");
					__SHOW(gen_box);
					
					//----------------------------------------------------------
					//-- get the #particles
					//----------------------------------------------------------
					const size_t gen_particles = size_t( Setup::GetReal(L,"gen_particles") );
					if( gen_particles <= 0 ) throw libc::exception( EINVAL, "gen_particles<=0" );
					__SHOW(gen_particles);
					
					//----------------------------------------------------------
					//-- get the temperature
					//----------------------------------------------------------
					const Real gen_temperature = Setup::GetReal(L,"gen_temperature");
					if( gen_temperature < 0 ) throw libc::exception( EINVAL, "gen_temperature<0");
					__SHOW(gen_temperature);
					
					//----------------------------------------------------------
					//-- create the frame
					//----------------------------------------------------------
					U3D cr;
					const size_t np = FCC::Count( cr, gen_particles, gen_box );
					frame.reset( new Frame( step->vectors(), np ) );
					
					//----------------------------------------------------------
					//-- place the particles
					//----------------------------------------------------------
					frame->set_box(gen_box);
					__SHOW(frame->particles);
					FCC::Build( (*frame)[Frame::POSITIONS], cr, frame->box, frame->ibox, 0 );
					
					//----------------------------------------------------------
					//-- assign masses
					//----------------------------------------------------------
					const Real gen_mass_min = Setup::GetReal(L,"gen_mass_min"); if( gen_mass_min <= 0 )          throw libc::exception( EINVAL, "gen_mass_min<=0");
					const Real gen_mass_max = Setup::GetReal(L,"gen_mass_max"); if( gen_mass_max < gen_mass_min) throw libc::exception( EINVAL, "gen_mass_max < gen_mass_min");
					__SHOW(gen_mass_min);
					__SHOW(gen_mass_max);
					const Real dm = gen_mass_max - gen_mass_min;
					for( size_t i=0; i < frame->particles; ++i )
					{
						const Real mass = gen_mass_min + Rand() * dm;
						(*frame)[ Frame::POSITIONS ][i].data = 1/mass;
						(*frame)[ Frame::VELOCITIES][i].data = mass;
					}
					
					//----------------------------------------------------------
					//-- generate velocities field
					//----------------------------------------------------------
					frame->initialize_velocities( gen_temperature );
					
					
					break;
				}
				
				
				throw exception("unknown frame_type=%s", frame_type.c_str());
			}
			assert( frame.__get()   != NULL );
			assert( pOutput.__get() != NULL );
			
			
			//==================================================================
			//
			// Prepare potential and dynamics
			//
			//==================================================================
			const Real epsilon = Setup::GetReal(L,"epsilon"); if( epsilon <= 0 )     throw libc::exception( EINVAL, "epsilon<=0");
			const Real sigma   = Setup::GetReal(L,"sigma");   if( sigma   <= 0 )     throw libc::exception( EINVAL, "sigma<=0"  );
			const Real cutoff  = Setup::GetReal(L,"cutoff");  
			__SHOW(epsilon);
			__SHOW(sigma);
			__SHOW(cutoff);
			
			if( cutoff > 0 )
			{
				//==============================================================
				//possible cutof
				//==============================================================
				if( cutoff  <= sigma ) throw libc::exception( EINVAL, "cutoff<=sigma");
				const Real cutout = Setup::GetReal(L,"cutout");
				if( cutout <= cutoff ) throw libc::exception( EINVAL, "cutout<=cutoff");
				
				//==============================================================
				// create the potential
				//==============================================================
				potential.reset( new LJCutOff(epsilon,sigma,cutoff) );
				
				//==============================================================
				// and the associated dynamics
				//==============================================================
				const string dynamics_type = Setup::GetString(L,"dynamics_type");
				__SHOW(dynamics_type);
				while( true )
				{
					//----------------------------------------------------------
					// full dynamics => discard the cutout...
					//----------------------------------------------------------
					if( dynamics_type == "full" )
					{
						std::cerr << "-- DISCARDING CUTOUT" << std::endl;
						dynamics.reset( new FullDynamics() );
						break;
					}
					
					//----------------------------------------------------------
					// single threaded simulation
					//----------------------------------------------------------
					if( dynamics_type == "cells:nothread" )
					{
						__SHOW(cutout);
						dynamics.reset( new CellDynamics( *frame, cutout ) );
						break;
					}
					
					//----------------------------------------------------------
					// multithreaded simulation, using YOCTO_THREADING
					//----------------------------------------------------------
					if( dynamics_type == "cells:threaded" )
					{
						__SHOW(cutout);
						dynamics.reset( new CellDynamicsMT( *frame, cutout ) );
						break;
					}
					
					throw exception("Unknowm dynamics_type=%s", dynamics_type.c_str());
				}
				
				
			}
			else 
			{
				//----------------------------------------------------------
				// single threaded simulation
				//----------------------------------------------------------
				std::cerr << "-- DISCARDING CELLS POSSIBILITY" << std::endl;
				potential.reset( new LJ(epsilon,sigma) );
				dynamics.reset( new FullDynamics() );
			}
			
			//==================================================================
			//
			// Get run time
			//
			//==================================================================
			const Real runtime = Setup::GetReal(L,"runtime");
			if( runtime < 0 ) throw libc::exception( EINVAL, "runtime<0");
			__SHOW(runtime);
			t_sim = runtime;
			
			//==================================================================
			//
			// Get kind of simulation
			//
			//==================================================================
			const string ensemble = Setup::GetString(L,"ensemble");
			while( true )
			{
				if( ensemble == "NVE" )
				{
					isothermal = false;
					break;
				}
				
				if( ensemble == "NVT" )
				{
					isothermal = true;
					break;
				}
				
				throw exception("Unknown ensemble %s", ensemble.c_str());
			}
			__SHOW(ensemble);
			
		}
		assert( potential.__get() != NULL );
		assert( dynamics.__get()  != NULL );
		
		
		
		////////////////////////////////////////////////////////////////////////
		//
		// Set up simulation
		//
		////////////////////////////////////////////////////////////////////////
		
		//======================================================================
		//
		// evaluate #ops per step
		//
		//======================================================================
		const size_t N          = frame->particles;
		const size_t num_eval_f = ( N*(N-1) ) >> 1;
		const size_t num_update = ( isothermal ? 10 : 4 ) * N;
		const size_t num_ops    = num_eval_f + num_update;
		
		//======================================================================
		// prepare progress indicator
		//======================================================================
		eta          prg;
		const size_t n_run = size_t( __ceil(t_sim)/(timings->dt_run) );
		Real         t_run = t_ini;
		
		prg.reset();
			
		//======================================================================
		// Starting forces
		//======================================================================
		dynamics->eval( *frame, *potential, isothermal );
		prg(1.0/n_run);
		
		std::cerr.flush();
		progress( prg, num_eval_f / prg.now());
		
		//======================================================================
		// initial output if necessary
		//======================================================================
		{
			frame->save( *pOutput, t_run, false );
		}
		
		////////////////////////////////////////////////////////////////////////
		//
		// Run Simulation
		//
		////////////////////////////////////////////////////////////////////////
		double ell = 0;
		for( size_t i=1; i <= n_run; ++i )
		{
			const bool   shall_save = 0 == (i%timings->every);
			const double ell_ini    = prg.now();
			step->predict(  *frame );
			dynamics->eval( *frame, *potential, isothermal );
			step->correct(  *frame, isothermal );
			ell += prg.now() - ell_ini;
			t_run = t_ini + timings->dt_run * i;
			
			if( shall_save )
			{
				const double nops = timings->every * (num_ops / ell);
				ell = 0;
				prg( double(i)/n_run );
				progress( prg, nops);
				{
					frame->save( *pOutput, t_run, true );
				}
			}
		}
		prg(1.0);
		progress( prg, 0);
		std::cerr << std::endl;
		////////////////////////////////////////////////////////////////////////
		//
		// Done...
		//
		////////////////////////////////////////////////////////////////////////
		
		
		
	}
	catch( const exception &e )
	{
		std::cerr << "*** " << e.what() << std::endl;
		std::cerr << "*** " << e.when() << std::endl;
		return -1;
	}
	catch(...)
	{
		std::cerr << "*** Unhandled exception!" << std::endl;
		return -1;
	}
	
	return 0;
}
